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Different formulations of Einstein’s equations used in numerical relativity can affect not only the 
stability but also the accuracy of numerical simulations. In the original Baumgarte-Shapiro-Shibata- 
Nakamura (BSSN) formulation, the loss of the angular momentum, J, is non-negligible in highly 
spinning single black hole evolutions. This loss also appears, usually right after the merger, in highly 
spinning binary black hole simulations, The loss of J may be attributed to some unclear numerical 
dissipation. Reducing unphysical dissipation is expected to result in more stable and accurate 
evolutions. In the previous work [I] we proposed several modifications which are able to prevent 
black hole evolutions from the unphysical dissipation, and the resulting simulations are more stable 
than in the traditional BSSN formulation. Specifically, these three modifications (Ml, M2, and M3) 
enhance the effects of stability, hyperbolicity, and dissipation of the formulation. We experiment 
further in this work with these modifications, and demonstrate that these modifications improve the 
accuracy and also effectively suppress the loss of J, particularly in the black hole simulations with 
the initially large ratio of J and the square of the ADM mass. 

PACS numbers: 04.25.Dm, 04.30.Db, 95.30.Sf, 97.60.Lf 


I. INTRODUCTION 

Development of numerical relativity has been rapid af¬ 
ter the breakthroughs in 2005 and 2006 (see, e.g., )- 

Numerical relativity has now become an indispensable 
and effective tool in the research of general relativity and 
relativistic astrophysics. It has been extensively studied 
in several areas; and applied to the construction of grav¬ 
itational waveform template banks for detectio n [5 11, to 
the kick phenomena of general binary systems [6j-|10], and 
to astrophysical problems such as the equation of state 
of neutron stars [lll4+| , electromagnetic counterparts of 
gravitational waves [15 . [l6| , gamma ray bursts [17] , ac¬ 
cretion disks |18] and so on. These applications and nu¬ 
merical investigations demand increasingly greater accu¬ 
racy, besides stability; thus it is important and necessary 
to unremittingly refine the formulations and schemes. 

Among the methods to enhance the stability and ac¬ 
curacy in numerical relativity, the 3 + 1 formulation of 
Einstein’s equations is favored by many researchers. The 
Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formu¬ 
lation [ll 0] is the most popular scheme, and it is usu¬ 
ally implemented in first-order-in-time and second-order- 
in-space finite-differencing codes. Many works have been 
focused on modifying the original BSSN formulation to 
achieve better numerical stability and accuracy f+fl f+l . 
For example, borrowing from ideas in the Generalized 
Harmonic (GH) formulation [2:, [27j], the Z4 conformal 
(Z4c) formulation [28| and the traceless-conformal and 
covariant Z4 (CCZ4) formulatio n [2fll| both show good 
constraint damping behavior 

There is still room to improve the BSSN formulation to 
obtain better stability and accuracy, e.g., see references 


[2lU23l . [3314351 ] . In [l[ , we adopt a different approach from 
the CCZ4 and Z4c formulations to modify the BSSN for¬ 
mulation. The basic idea is to suppress the numerical 
error by adding combinations of constraint terms to the 
field equations in the BSSN formulation without chang¬ 
ing the solution analytically, to modify the leading terms 
of the field equations. And we demonstrated that our 
modifications achieved more stable simulations than the 
traditional BSSN formulation. Specifically, our last work 
fl] has shown improvements in constraint damping and 
in the late-time behavior of the gravitational waveforms. 
In this work, we would like to emphasize the effectiveness 
of these modifications on the evolution of the black holes 
with higher spins, hoping to meet the demand of model¬ 
ing extreme sources for the gravitational wave detection. 

It was found in [36] that the angular momentum de¬ 
cays right after the final black hole forms in the binary 
black hole evolution simulations. For the single spin¬ 
ning black hole, the angular momentum also decays when 
the dimensionless spin s/m 2 > 0.75 (compare the cases 
s/m 2 = 0.53 and s/m 2 = 0.9 in Fig. 4 of [36|). This 
decay is neither due to the resolution nor the initial sep¬ 
aration of the binary 0]. In contrast, the result from the 
SpEC code does not show this tendency [37]- For com¬ 
paring those results, we plot in Fig. [T] x/ = s//m 2 of the 
final black hole as a function of the initial spin parame¬ 
ter Xi = s i/ m 'i for individual black hole component. We 
find that they are consistent when Xi < 0-75 and, when 
Xi getting larger, Xf i n the traditional BSSN formula¬ 
tion becomes smaller than that in the GH formulation. 
The difference should not be attributed to whether the 
spectral method or the finite-differencing method is used. 
However, it is yet unclear if the issue comes from the for- 
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FIG. 1: Dimensionless spin parameter \f °f the final black 
hole as a function of the initial dimensionless spin parameter 
Xi of the individual black hole in the binary black hole evo¬ 
lution for the BSSN and GH formulations. The data ’BSSN’ 
given in [36j] are calculated through the BAM code with BSSN 
formulation and the finite-differencing method. The data 
’GH’ given in [37|] are calculated through the SpEC code with 
GH formulation and the spectral method. The data ‘BSSN 
before decay’ given in [36[ is calculated through the LEAN 
code [3(f]. It corresponds to the result right after the final 
black hole forms, and the data ‘MODBSSN’ are calculated 
through AMSS-NCKU code in current work, and the modi¬ 
fied BSSN formulation (check the main text for more detail) 
is used. 

mulation itself or from the puncture method. In this 
work, we therefore try to resolve this problem by simu¬ 
lating single and binary black hole evolutions with the 
modified BSSN formulation as proposed in |T}. We will 
show that the angular momentum is more accurate and 
its conservation is much better than in the traditional 
BSSN formulation. Via the better conservation of the 
angular momentum, it is also expected that the accuracy 
of the other related physical quantities will be improved 
at the same time with the modified formulation. 

The rest of this work is organized as follows: In the 
next section, we will give an explicit description of the 
modifications to the BSSN formulation, discuss the re¬ 
lated accuracy problems, and describe the numerical im¬ 
plementation. We then report on the test results on sin¬ 
gle spinning black hole in Sec. lIII Al The results for highly 
spinning binary black hole are presented in Sec. IIII Bl 
And the discussion and summary will be presented in 
the Sec. UV] Throughout the paper, geometric units with 
G = c = 1 are used. Einstein summation rule is adopted 
unless stated explicitly. 


II. MODIFICATIONS AND NUMERICAL 
IMPLEMENTATIONS 

The BSSN formulation and the numerical recipes for 
implementation have been described in details in previ¬ 


ous articles (HI, H3- Here we only mention several major 
steps that have usually been adopted [21, [33j in the tra¬ 
ditional BSSN formulation: 

• In order to enforce the algebraic constraints of the 
unimodular determinant of 7^, i.e., 7=1, and 
of the tracelessness of Aij, i.e., 7 ^Aij = 0, the 
numerical values of 7y and Aj 3 are replaced with 
7 ij -A 7~ 1/3 7ij, Ai 3 -a A^ij) after every time step, 
wherein the two indices in the angle bracket () is 
taken to be its symmetric and traceless part. 

• The conformal connection functions T* are pro¬ 
moted to be independent variables in the BSSN 
formulation, which leads to the T-constraints Q % = 
f’ —fg = 0, where fg = 7 jfc f l jfe . The conventional 
approach to enforce the T-constraints is to replace 
all the undifferentiated T* with fg. 

• The high-order Kreiss-Oliger (KO) method is em¬ 
ployed to dissipate effectively the numerical noise. 

These traditional approaches with suitable gauge condi¬ 
tion have enabled fruitful studies on the black hole prob¬ 
lem. Yet earlier investigations, e.g., [36], indicated that, 
in some near extreme situation, the traditional BSSN for¬ 
mulation is not robust enough to conserve the constraints 
and global quantities. We plan to test the following modi¬ 
fications which have been introduced in |lj, and compare 
the results from our modifications with those from the 
traditional BSSN formulation. The three proposed mod¬ 
ifications are as follows: 


A. Modification Ml 

Instead of replacing all the undifferentiated f l with fg, 
Ml modifies the conformal 3-connection appearing in the 
right-hand-side of all the field equations, and changes 
the linear terms in the field equation of TL The new 
conformal 3-connection in all field equations now takes 
the form 

fv - (1) 

where T) = T k k i = (ln-y/7)^ vanishes analytically, but 
could be nonzero due to numerical error. This expression 
is motivated by the unique algebraic decomposition for 
any third-rank symmetric tensor with two indices. See 
[l] for the details. 

To change the behavior of the linear term in the field 
equation of f% we replace the original field equations of 
f 1 with 

d t r =2a[r jk A^ - + 6 A»(t>j] - 2 A^ a 3 

+ ptr d - f '.T y + f k p\ jk + \i ij P k ,j k 
+ |(/3 fc jfc - 2 aK)r - (1 + £)Q(Y)Xg\ (2) 
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wherein ©(a;) is the step function 


0(:r) 


0 if x < 0 
1 if x > 0 


and X 1 is 



2aK) - 



(3) 

(4) 


Note that the index with hat, i.e., i. means that no index 
summation is carried out with respect to this index. £ is 
chosen to be 1 in all cases in this work. This modifica¬ 
tion plays an indispensable role in the whole modification 
scheme to enhance both the stability and accuracy of the 
system. 

Notice that there is one term in eqn 0 including a 
step function, i.e., (1 + ftQ(X l )X l G l . Due to its switch 
character and the possible sign fluctuation of the numer¬ 
ical value of its argument A 1 when A* is close to zero, the 
step function should be sensitive to the resolution used 
in simulations. So we expect that this modification could 
affect majorly the numerical convergence behavior of the 
modified BSSN formulation. 


B. Modification M2 

The idea behind M2 is similar to that in obtaining 
Eq. (U). The algebraic structure of dt'jij, similar to the 
algebraic structure of the conformal 3-connection, allows 
us to write the 7y-field equation as 


dftij t fA ~(ij T u 3 (, Qj j ^ Xij b Qk) (a) 

and we set a = 1/10 in this work. This modification 
enhances the hyperbolicity of the system and propagates 
the constraint violation residual away effectively. 


C. Modification M3 


m might change the convergence order of a system to 
be only second-order accurate at most since the addtion 
term in the equation, i.e., is only proportional 

explicitly to h 2 . However, this is not the case. If one 
system is p th -order convergent, then the momentum con¬ 
straint A4 i ~ 0 will converge to zero with the rate of h p . 
So will the term M uj\. Therefore, if we combine the 
convergence order of and the the multiplier A 2 , 

the term introduced in eqn ([6]) will converge to zero with 
the rate h p+2 , which convergence is faster than the rest 
of the system. Thus this modification will not reduce the 
convergence order of the system analytically. 


D. Numerical Implementation 

The AMSS-NCKU code with the standard moving box 
style mesh refinement [H H3, [H| is used in this work. We 
used 10 mesh levels, all of which are fixed in the cases 
of single black hole evolution, and the finest 3 levels are 
movable in evolving the binary black holes (BBHs). In 
each fixed level, we used one box with 128 x 128 x 64 grids 
with assumed equatorial symmetry. The outermost phys¬ 
ical boundary is 512M and this makes the finest resolu¬ 
tion to be h = M/64. For the movable levels, two boxes 
with 64 x 64 x 32 grids are used to cover each black hole. 
In time direction, the Berger-Oliger numerical scheme is 
adopted for the levels higher than four. 

The moving puncture gauge condition 

d t a = fta^ - 2 aK, (7) 

d t ft = ^B i +ftft J , (8) 

d t B l = dft 1 - r]B l + /3 j Bj - ft fV. (9) 

is used and has been shown to give good behavior for 
the black hole simulations in [34[. In this paper we use 
r] = 2 M with M being the ADM mass of the given con¬ 
figuration. 


This dissipation type of modification M3 is motivated 
from [23|. The major difference is that we use the sym¬ 
metric traceless part of the partial derivative of the mo¬ 
mentum constraint (instead of the symmetric part of its 
covariant derivative as in [23[) to re-write the A^-field 
equation as 

dtftj —► dtftj + (6) 

wherein A it is the momentum constraint, h is the grid 
width. This modification provides a dissipation mecha¬ 
nism on Aij, and serves as a natural alternative to the 
KO dissipation. In this work, we apply this modifica¬ 
tion instead of the KO method to check its capability in 
dissipation and also compare its effect with KO’s. 

There is a concern about the convergence of the whole 
system with this modification. At first glance, Equation 


III. NUMERICAL RESULTS 
A. Single black hole tests 

In this subsection, we test our modifications in spin¬ 
ning single black hole (SBH) cases with respective ini¬ 
tial dimensionless spin parameters xo = Jq/AIq = 0.53, 
0.9, and 0.923. To generate these sets of puncture initial 
data with unit ADM mass, the bare mass and the spin 
parameter in the 2-direction are set to be m = 0.872335, 
0.3528, and 0.215898 and = 0.53, 0.45, and 0.472466 
respectively, as the input for the TwoPuncture solver. 
Note that xo = 0.923 is nearly the maximal spin that 
the conformally flat Bowen-York initial data can achieve 
[39| . The global quantities such as the ADM mass M and 
the angular momentum J are calculated with the surface 
integrals at R = 50M, as described in [34j]. 
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FIG. 2: Dimensionless spin parameter x as a function of time for xo = 0.53 (left column) and xo = 0.9 (right column) in the 
single black hole evolutions. The modified BSSN formulation (solid red line, marked as MODBSSN) is shown to produce less 
noise in the spurious radiation as well as preserving x better than the traditional BSSN formulation (dashed line, marked as 
BSSN) in the higher spin case. 


Single BH, = 0.923 Single BH, % 0 = 0.923 



FIG. 3: Left: Dimensionless spin parameter \ as a function of time for xo = 0.923 in the single black hole evolution. The 
modified BSSN formulation (solid red line, marked as MODBSSN) is shown to preserve x better than the traditional BSSN 
formulation (dashed line, marked as BSSN) in this near extreme case. Right: Power spectrum of the corresponding data plotted 
in the left panel. 


The dimensionless spin parameter x = J/M 2 for the 
single black hole simulation is shown in Fig.[2]and Fig. [3] 
For the cases with the traditional BSSN formulation, it 
shows that the result is consistent with the Fig. 4 of 1361 ]. 
It is also clear in the figures that the proposed modifica¬ 
tions greatly reduce the overall noise level, and diminish 
the fluctuation before t = 200 for the high-spin case. The 
curves of \ with the modified BSSN formulation show less 
decay in each case. For the lower spin case, xo = 0.53, 
the curve for the modified BSSN formulation (red solid 
lines) is basically the same as the one for the traditional 
BSSN formulation (dashed lines). As the spin becomes 
higher, the loss of the angular momentum is more severe. 
For the xo = 0.9 case, the dimensionless spin drops more 
than 11% to 0.7986 at t = 1200 in the traditional BSSN 
formulation, compared to the modified one in which x 


drops less than 5%. This result indicates that the modi¬ 
fied BSSN formulation is capable to conserve the angular 
momentum better than the traditional BSSN, especially 
for the high spin cases, i.e., xo >0.75. 

The time evolution of the dimensionless spin is shown 
in the left panel of Fig.[3]for the near-extreme single black 
hole with xo = 0.923, which is nearly the maximal value 
that the conformally flat Bowen-York data can achieve. 
The dimensionless spin drops about 15% to 0.79 at t = 
1000 in the traditional BSSN formulation, compared to 
the modified case in which the change of x is less than 
5%. It shows that the modified BSSN formulation is more 
effective in conserving the angular momentum over the 
traditional BSSN formulation, even in the fast-spinning 
SBH case. 

It is interesting to study the different effect on the sim- 
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FIG. 4: The puncture trajectory of the binary black hole evolution with initial — 0 (left) and \i = 0.9 (right). The com¬ 
parisons between the result with the traditional BSSN formulation (dashed line) and that with the modified BSSN formulation 
(solid red line) are shown. 




FIG. 5: The £ = 2, m = 2 mode of the Newman-Penrose scalar calculated at R = 50M, with initial = 0(left) and 
Xi = 0.9(right). This mode is the major component of gravitational radiation during the merger in the binary black hole case 
with spin parallel to the orbital angular momentum. Only the real part of 4G is shown here. These two waveforms are almost 
identical for the traditional BSSN formulation (dashed line, marked as BSSN) and the modified one(solid red line, marked as 
MODBSSN). The binary composed of higher-spin hole spend more inspiral period due to the spin hang-up effect. 


illations between the 5 th -order KO dissipation and Mod¬ 
ification M3. From a naive observation on Fig. [2] and 
the left panel of Fig. [3J we found that the KO method is 
good at eliminating relatively higher frequency numeri¬ 
cal noise. It can be seen that the result in right panel of 
Fig. [2] for the xo = 0.9 case with the traditional BSSN 
formulation (dashed line) is smoother than its counter¬ 
part with the modified BSSN formulation (red solid line), 
despite the spin drop in the former one. In contrary, the 
lower frequency numerical noise appearing in the tradi¬ 
tional BSSN formulation is diminished significantly with 
the modified BSSN formulation. This can also be seen in 
the left panel of Fig. [5] wherein the two major fluctuations 
at t pa 580 and t ~ 1150 with the traditional BSSN formu¬ 
lation (dashed lines) disappear with the modified BSSN 


formulation. It also can be seen in the right panel wherein 
the severe fluctuations before t ss 150 with the traditional 
BSSN formulation is effectively suppressed with the mod¬ 
ified BSSN formulation. It has similar behavior in the left 
panel of Fig. [3] 

To understand this phenomenon better, a Fourier anal¬ 
ysis method is applied to the xo = 0.923 single black hole 
case. In the right panel of Fig. [3] we show the correspond¬ 
ing power spectrum of the data in the left panel. From 
this power spectrum, we can see that the KO method in 
the traditioinal BSSN formulation only dissipates some 
high-frequency (/ ~ 0.34 — 0.43) noise better. For the 
noise in the most other frequency, the M3 method in the 
modified BSSN formulation is much more efficient in dis¬ 
sipation. This difference results in the different behavior 
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Black hole binary, X; = 0 Black hole binary, X; = 0.9 




FIG. 6: ADM mass M, angular momentum J, and dimensionless spin parameter \ as functions of time for the initial = 0 
(left column) and the \i = 0-9 (right column) in the binary black hole evolutions. The modified BSSN formulation (solid red 
line, marked as MODBSSN) has better control than the traditional BSSN(dashed line) on the early noise level and the spin 
drop after merger for the higher spin case. 


in the left panel as we can see. The line for traditional 
BSSN formulation has larger amplitude oscillations with 
the intermediate frequency. The similar results can be 
seen in Fig. [2] The above result indicates that the KO 
dissipation and M3 suppress the numerical noise in differ¬ 
ent frequency ranges. It is noted that usage of M3 does 
not introduce any artificial dissipation and thus the field 
equation of Ay is analytically equivalent to the original. 


B. Binary black hole tests 

In this subsection, we apply our modifications to the 
equal-mass black hole binary. Each black hole in the 
binary has the spin aligned with the orbital angular mo¬ 
mentum and the dimensionless spin parameter Xi — 0-9 
initially. As the reference, we also run an equal-mass 
BBH with Xi = 0 for each black hole. The initial pa¬ 
rameters for each hole as the input of the TwoPuncture 
solver are listed in Table Q] 

TABLE I: Parameters for the binary black hole puncture ini¬ 
tial data 


Xi 

0 

0.9 

bare mass 0.483 

0.1764 

f 

±3.257y 

±2.966y 

P 

=F0.133x q=0.12616x 

s 

0 

0.225z 


Firstly we would like to check if our modifications give 
any changes in these well-tested BBH cases, compared 


to the traditional BSSN formulation. Figures |4] and [5] 
show the almost identical puncture trajectories and the 
£ = 2, m = 2 mode of the Newman-Penrose scalar 'F4 
at R = 50 M for the cases of Xi = 0 and Xi = 0-9. The 
£ = 2, m = 2 mode gives the major component of grav¬ 
itational radiation during the merger in the case of the 
binary black hole with spin parallel to the orbital angular 
momentum. Here we only show the real part of ^4. This 
result is expected in developing new modifications since 
all of these formulations are analytically equivalent to 
Einstein’s field equations and should give same physics. 

The ADM mass M, the angular momentum J, and the 
dimensionless spin parameter x are shown in Fig. [6] for 
the BBH cases with the initial dimensionless spin \i = 0 
(left) and Xi = 0.9 (right). After the merger at t = 250 in 
the Xi = 0 BBH case, M and J decrease by 4% and 27% 
respectively due to the gravitational radiation. And the 
dimensionless spin parameter after the merger is x = 0.68 
at t = 250 to x = 0.67 at t > 800. Thus x decreases less 
than 2% after t = 250 until the end of simulation. It 
also shows that the modified and traditional BSSN for¬ 
mulations give the same result (in the left panel) in the 
initially slowly spinning BBH case. In the Xi = 0.9 case, 
after the merger at t = 350, the gravitational radiation 
decreases the value of M and J by 8% and 39% respec¬ 
tively. As shown in the right panel for the Xi = 0.9 case, x 
in the traditional BSSN formulation decreases more than 
10% from t = 350 to 800 (% w 0.75 in our extended run 
for t > 1900). This result from the traditional BSSN for¬ 
mulation is consistent with the discovery in [36j] in which 
the final J will decay considerably for Xi > 0.75. For the 
Xi = 0.9 case with the modified BSSN formulation, x de- 
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creases only by 1% from t = 350 to 800. The decrease is 
still less than 2% in the extended run for t > 1900. The 
results in the BBH cases, combined with that in the SBH 
cases, indicate clearly that our modifications handle the 
highly spinning black holes much better than the tradi¬ 
tional BSSN formulation, while yielding the same results 
as in traditional BSSN formulation in the slow spinning 
black hole cases. 


C. Numerical Convergence 



FIG. 7: Effect of the extraction radius on the ADM mass 
integration. The plot corresponds to the spinless binary black 
hole case. The r’s in the legend are the extraction radii used 
in this case. 

For both the gravitational wave extraction and the cal¬ 
culation of the global quantities, they are numerically 
integrated on the sphere of radius r = 50M. This finite 
radius for integration could affect the accuracy of the am¬ 
plitude of the gravitational waveform \I/4. However, the 
effect from the integration sphere of finite radius should 
be roughly the same with either the traditional BSSN 
formulation or the modified one on any case. Since we 
are only concerned about the relative difference between 
these two formulations, the effect from the extraction ra¬ 
dius becomes unimportant. For the global quantities, 
e.g., the ADM mass and the angular momentum, the 
integration sphere of finite radius, e.g., r = 50 M, may 
result in some artificial drift as shown in Fig. [G] Never¬ 
theless, when an integration sphere with larger radius is 
applied to the case, such kind of drift diminishes. Here 
we use the spinless binary black hole case as an exam¬ 
ple in Fig. [7] for illustration. According to Fig. [7] the 
result with r = 50M are basically same as the ones with 
r = 80 M and r = 120M during the merger phase. The 
drift in the case with r = 50M only shows during the 
ringdown stage. And the drift can be easily diminished 
with larger radii, e.g., r = 120M in Fig. 0 However, in 
order to compare our result with the one in [36] closely, 
in this work we still take the radius r = 50M, same as 
used in [36j . 

The physical boundary used in simulations may also 
affect the gravitational wave and the calculation of the 
global quantities. In order to investigate such possi¬ 


ble effects, we have tested the simulations with farther 
boundaries. And our results show that the effect from 
the boundary condition is ignorable in the current work. 

According to the arguments in Sec. [111 analytically we 
do not expect Modification M3 to affect the convergence 
of the system, and we do expect that Modification Ml 
definitely affects the system’s numerical convergent be¬ 
havior due to the switch character of the step function 
in Eq. ©. Here we would like to both check the con¬ 
vergence order of the system with the modified BSSN 
formulation and verify these arguments numerically. 

Firstly, we show the system’s convergence with our 
modifications in the left-side panels of Fig. [51 compared 
with the one with the traditional BSSN formulation. 
Here we use the Xi = 0 binary black hole case as an exam¬ 
ple for demonstration. In the two plots we study the con¬ 
vergence of the phase and the amplitude for the gravita¬ 
tional wave respectively. According to the plots, the tra¬ 
ditional BSSN formulation results in overall 3.3 th -order 
convergence for the system in both the phase and am¬ 
plitude of ^ 4 , which is roughly consistent with the ideal 
convergence of fourth-order with the numerical method 
used in this work. On the other hand, it shows in the 
panels that the modified one results in only first-order 
convergence for the system. Since we already expect that 
some of our modifications will affect the convergence of 
the system, the result is understandable, although the 
order of convergence is still considered low. And we can 
see from Fig. [5] that the numerical error with the modi¬ 
fied BSSN formulation is much smaller than the one with 
the traditional BSSN formulation, especially in the lower 
resolutions. This merit for the modified BSSN formu¬ 
lation could compensate for its disadvantage of having 
lower order convergence. And the convergence behavior 
showed in Fig. [5] is general for all the cases we have done 
in this work. 

Secondly, we would like to confirm the theoretical anal¬ 
ysis that it is Modification Ml, not M3, in the modified 
BSSN formulation which majorly affects the convergence 
order. By using again the \i = 0 binary black hole case 
as an example, we show the result in the right-side panels 
of Fig. [8] When we apply the traditional BSSN formu¬ 
lation + Ml to the case, the resulted convergence order 
is first-order, which is roughly the same as the conver¬ 
gence order for the case with the whole modified BSSN 
formulation. Meanwhile, when we apply the traditional 
BSSN formulation + M3 to the case, the resulted con¬ 
vergence order is 2.5 th -order, which is a little lower than 
the case with the pure traditional BSSN formulation, but 
quite higher than the one with the modified one.. The 
result tells that Ml is the key modification which affects 
majorly the convergence behavior of the system, as we 
expected. However, it also shows that M3 lowers minorly 
the convergence order of the system. This indicates that 
the convergence order of M3 might not be numerically 
as g ood as the expectation from our analytical argument 
[40( • In conclusion, the cases with our modifications have 
first-order convergence, which is lower than the one with 
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FIG. 8 : Convergence of gravitational wave for the \i = 0 binary black hole case. The left-side panels show the phase differences 
and the amplitude differences of 'I '4 respectively between the high and medium resolutions (solid line), and the medium and low 
resolutions (dotted line), with both the traditional BSSN formulation (BSSN, black line) and the modified one (MODBSSN, 
red line). Here we use f(p) to denote the factor of order p. The right-side panels show the phase differences and the amplitude 
differences of ' 3/4 respectively between the high and medium resolutions (solid line), and the medium and low resolutions (dotted 
line), with the traditional BSSN formulation+Modification M3 only (M3, black line) and the traditional one+Modification Ml 
only (Ml, red line). 


the traditional BSSN formulation, but our modifications 
give more accurate result than the traditional BSSN for¬ 
mulation at a given resolution. And Modification Ml is 
the key factor in our modifications to affect the conver¬ 
gence behavior. 


IV. DISCUSSION AND SUMMARY 

In this work, we applied our modifications of the BSSN 
formulation to study the total angular momentum con¬ 
servation issue in black hole evolutions with the standard 
Bowen-York puncture initial data. We found that the 
non-negligible loss of angular momentum for highly spin¬ 
ning black holes mentioned in [36| can be greatly cured 
with our modifications. The improvements are obvious 
for near-extreme cases, as in the SBH case shown in Fig. [2] 
and Fig. [3] and the BBH case in Fig. 0 It has also been 
shown in the previous section that the modified BSSN 
formulation does not introduce any unphysical effects. 
Improving the conservation of the angular momentum 
usually leads to certain improvement on the accuracy of 
the results in black hole evolutions. Therefore we expect 
that our modifications will provide better performance 
in black hole evolution simulations than the traditional 
BSSN scheme. 

Modification Ml is the most important to the conser¬ 
vation of the angular momentum since the field equation 


of the conformal connection function f * is closely related 
to the (angular) momentum vector. We find that Eq. © 
and setting £ = 1 gives quite robust and stable runs. 
But due to the switch character of Ml and the possible 
sign fluctuation of the numerical value of its argument 
when the argument is close to zero, the step function is 
sensitive to the resolution used in simulations. So this 
modification affects majorly the numerical convergence 
order of the modified BSSN formulation. M2 is able to 
enhance the hyperbolicity of the system, especially for 
the evolution of Tb However, its mechanism and the 
optimal choice of a need further investigations. 

Instead of the KO dissipation method used in the tra¬ 
ditional BSSN formulation, Modification M3 is used in 
the modified BSSN formulation in this work. We can see 
in Sec. lIIIl that M3 is able to diminish effectively some in¬ 
termediate frequency noise with larger amplitude; while 
the KO dissipation is good at eliminating the higher fre¬ 
quency numerical noise. To some extent, Modification 
M3 is complementary to the KO method in dissipating 
numerical error. However, the advantage of M3 is that it 
comes from the derivative of the momentum constraint. 
Thus, applying M3 to the BSSN formulation is always 
legitimate and safe as long as the momentum constraint 
holds. In contrast, the application of the KO method 
is not always safe since it is an artificial addition to the 
field equation, although it is convenient and effective in 
numerical relativity. It is possible that usage of the KO 
method leads to deviations of the numerical result from 
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the solution hypersurface, especially when the result is 
sensitive to the initial data. It will be a good idea to use 
both Modification M3 and the KO method in dissipating 
the numerical error in the future simulations. 

As mentioned at the beginning of Sec. [TTJ in order to 
enforce the algebraic constraints of the unimodular deter¬ 
minant of 7 ij, and of Ajj being traceless, the numerical 
values of 7^ and A^ are replaced with 7 —> 7^ 1//3 7ij, 
Aij —> A(ij\ after every time step in the traditional BSSN 
formulation. However, the modification 


Izz 


A 


yy 


1 + lyyllz ~ ^Ixylyzlxz + lxx% z 
Ixxliyy ~ 7 xy 

K* + A z z + A xy y*v + AyrfV* 

■y yy ’ 


( 10 ) 

( 11 ) 


is employed instead in jT] to enforce the two constraints. 
The results in [l] have shown that this modification gains 
better stability compared to the traditional BSSN formu¬ 
lation. Differing with [l]], in this work we use the tradi¬ 
tional recipe of 7 ij —> 7 _1 ^ 3 4 5 6 7 8 9 10 7ij, A ij -A Auj\ instead of 
applying the modification Eqs. (Hhd & dm). This is be¬ 
cause the denominators in Eqs. (TtOlTTl) can be very small 
near the singularities and thus the numerical values of 
the replaced 7 zz and A yy can have unexpected fluctua¬ 
tion which can crash the code. However, this modifica¬ 
tion can still be applied to the BSSN formulation if there 
is no singularity or if an excision method is used in the 
black hole evolution simulations. 

In the modifications, we introduce some terms related 
to the spacial resolution used in the numerical simulation. 
These terms reduce the convergence order from fourth- 
order to first-order in our implementation. But compared 


with the traditional BSSN formulation, the numerical er¬ 
ror resulted in the modified BSSN formulation is still ob¬ 
viously smaller than the one from the traditional BSSN 
formulation with a reasonably fine resolution. 

In this work, we demonstrated the simulations with 
the modified BSSN formulation in which the angular mo¬ 
mentum conservation is better than in the traditional 
BSSN formulation. Thus this modified BSSN formula¬ 
tion should improve the accuracy in the punctured black 
hole evolutions. Our modifications are imposed on the 
field equations of the physical variables 7 ij, A y and T 1 2 , 
instead of the gauge variables a and (3 l . Therefore, we ex¬ 
pect that the modified BSSN formulation can be applied 
generally to various scenarios to give improved results in 
numerical relativity. 
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